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We link linear prediction of Chebyshev and Fourier expansions to analytic continuation. We push 
the resolution in the Chebyshev-based computation of T = 0 many-body spectral functions to a 
much higher precision by deriving a modified Chebyshev series expansion that allows to reduce 
the expansion order by a factor ~ We show that in a certain limit the Chebyshev technique 
becomes equivalent to computing spectral functions via time evolution and subsequent Fourier 
transform. This introduces a novel recursive time evolution algorithm that instead of the group 
operator e~ lHt only involves the action of the generator H. For quantum impurity problems, we 
introduce an adapted discretization scheme for the bath spectral function. We discuss the relevance 
of these results for matrix product state (MPS) based DMRG-type algorithms, and their use within 
dynamical mean-field theory (DMFT). We present strong evidence that the Chebyshev recursion 
extracts less spectral information from H than time evolution algorithms when fixing a given amount 
of created entanglement. 


I. INTRODUCTION 

Expanding the spectral density A(co) of an operator H 
in the monomes cu n via the moments 

Pn° n = J du>A(u>)u n , 

is a tool that originates in the early days of quan¬ 
tum mechanics. Computing these moments iteratively 
though is numerically unstable ’ and one replaced ex¬ 
pansions in uj n by expansions in such polynomials p n (co) 
of degree n that can be stably computed. ’ A prominent 
example for p ra (w) are Chebyshev polynomials, whose as¬ 
sociated three-term recursion is stable as it does not ad¬ 
mit a so-called minimal 1 solution. 

After the development of stable recursions the next 
step in the mid 1990s was the introduction of kernels that 
damp the erroneous Gibbs oscillations of truncated poly¬ 
nomial expansions of discontinuous functions,' “ which 
lead to the kernel polynomial approximation. It deals 
with redefined series expansions that represent the con¬ 
volution of the expanded function with a broadening ker¬ 
nel, like a Gaussian or Lorentzian. This technique has 
been reviewed in Ref. and more recently in Ref. 9 from 
a numerical linear algebra perspective. 

In this paper, we drop the idea of such broadening ker¬ 
nels in frequency space or the equivalent damping or win¬ 
dowing kernels in the associated Fourier or Chebyshev 
expansions. Instead, we employ the fundamentally dif¬ 
ferent technique of linear prediction. Linear prediction 
is a linear recursive reformulation (Appendix C) of the 
non-linear problem to fit the surrogate function 

g(t) = y aje tu>it , a„w,eC, lei, (1) 

i 

to given numerical data {t n ,g n }. Due to linearity, linear 
prediction is able to treat superpositions of hundreds of 


terms, and by that reliably extracts much information 
about an underlying function from its local knowledge 
{tmgn}- In order for this to be meaningful, the underly¬ 
ing function, e.g. a Green’s function, must be compatible 
with (1). 

In particular, we note that (1) can serve as an ansatz 
for analytic continuation of a zero-temperature Green’s 
function 

G(t) = (2) 

where \ifo) is a single-particle excitation of the ground- 
state | Eq) of H, for example the creation of a fermion 
j^o) = c'\Eq). Note that in the case of fermions, 
(2) describes only the t > 0 contribution (then usu¬ 
ally more precisely denoted G > (t)) of the full fermionic 
Green’s function. G(t) is analytic everywhere in the com¬ 
plex plane except for t — i ioo and thereby allows for 
an analytic continuation of G(t) from a local descrip¬ 
tion {t n ,G(t n )} to the domain [f 0 ,oo). This analytic 
continuation is highly different from the ill-conditioned 
problem of continuing the frequency-space represented 
Green’s function from a domain in the complex plane 
(e.g. the imaginary-frequency axis or a parallel of the 
real-frequency axis) to the real-frequency axis, where the 
frequency-space Green’s function has poles. 

In the context of Green’s functions, linear prediction 
has for the first time been used to extrapolate the time 
evolution of the spin structure factor in the one dimen¬ 
sional Heisenberg model. ’ While for the spin-1 model 
it was clear that the ansatz (1) is justified as the time 
evolution is dominated by a small number of magnons 
whose excitation energies correspond directly to the fre¬ 
quencies u>i in (1) , 1 this was not the case for the spin-1 
model. In the latter, spinons dominate which lead to 
an (infinitely) high number of poles on the real-frequency 
axis, and the direct correspondence of pole energies and 
frequencies Ui in (1) is lost. Still the ansatz works' in 
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an approximate sense by extracting effective frequencies. 

For the computation of spectral functions, the use of 
linear prediction for the time evolution of Green’s func¬ 
tions provides a highly attractive alternative approach to 
the usual damping or windowing in real-time or broad¬ 
ening in frequency space: an approach that enhances 
resolution in frequency space. Up to now, it is not en¬ 
tirely clear in which cases this is controlled. On the other 
hand, the approach of damping the truncated series ex¬ 
pansion cannot be considered controlled, too: Although a 
broadened function /^(w), which is for Gaussian broaden¬ 
ing given by /» = ^ / du'e~^'^ con¬ 

verges uniformly to the underlying original function f(co) 
for 77 —> 0 , extraction of information ( deconvolution ) from 
about /(w) is uncontrolled as it corresponds to the 
problem of analytic continuation from a domain in the 
complex plane to the real axis. 

Recently, Ref. 13 suggested to extrapolate the Cheby- 
shev expansion of a spectral function using linear predic¬ 
tion, albeit only justified by the empirical success. In the 
remainder of this introduction, we place these results in 
the context of the preceding discussion, and by that put 
this approach on more firm grounds. 


A. Chebyshev and Fourier transformation basics 


The Chebyshev polynomials of the first kind 

T n {x) = cos (narccos(:r)) (3) 

can be generated by the recursion 

T n (x) = 2xT n _i(x) - T n _ 2 (x), Xj = x, T 0 = 1, (4) 

which is numerically stable if |a;| < 1. Chebyshev polyno¬ 
mials are orthonormal with respect to the weighted inner 
product 


dx w n (x)T m (x)T n (x) = 5 nm , 


.w = 


2 - 6 , 


nO 


r\/l — . 


(5a) 

(5b) 


Any integrable function /(aOUef-i . 1 ] can be expanded in 
T n (x) 


Analogously, any integrable function /(w)| wg [_| j |], 
where a £ R, can be expanded in a Fourier series 

OO 

/M = ^ E e” 4 "/(tn), (7a) 

n=— 00 
/•“/2 

f{tn)= / dwe _lwt "/(w), t n = ~, (7b) 

J —a/2 a 

which represents a Fourier transform for a —» 00 . 


B. Expansion of a spectral function 

Consider now the expansion of the spectral function 
A(u>) of a Hamilton operator H with respect to a refer¬ 
ence energy E le f and a state |^o) as i n (2) 

A(w) = (il’o\S (uj - (H - E ie f))\ijjo). (8) 

The spectral function is related to the Green’s function of 
(2) via its Fourier transform: A{lS) = -UmG(w + i0 + ). 

The coefficients of the Fourier expansion can be com¬ 
puted by inserting an identity of eigenstates JA \Ei)(Ei\ 
in the integral over the delta function S(oj — (H — E re f)) 

f(t n )=[ dxve~ %utn A(u) = (ip 0 \ip(t n )), (9a) 

J-a/2 

\i>(t n )) = e~^ H ~ ET ‘ t)tn \ V> 0 ), t n = (9b) 

a 

In order for (9) to hold true, a must be chosen large 
enough for that the support of A(uf) is contained in 
[—1,|]. A sufficient condition for that is spec(7f — 
E re {) C [— |, |], which is possible as we consider opera¬ 
tors H with bounded spectra. Eq. (9b) makes it obvious 
that a has the meaning of an inverse time step. 

To compute the coefficients for the Chebyshev expan¬ 
sion, we need to consider a spectral function whose sup¬ 
port is contained in [—1,1], For this, introduce a rescaled 
and shifted version of H with appropriately chosen con¬ 
stants a and b 

'H a ,b = — —— +b , x=-+b, (10) 

a a 

where a can again be considered an “inverse time step” 
and T~L is dimensionless. Note that in Ref. 14, the defini¬ 
tion of b differed from the one here by a factor a. Then 

A a ,b(x) = {%\ 6 (x - 74 , 6 ) 1 ^ 0 } ( 11 ) 


OO 

f(x) = y^ j w n (x)p, n T n (x), ( 6 a) 

n—0 

Pn = J dxT n (x)f(x ), ( 6 b) 

where the definition of the so-called Chebyshev moments 
p n via the non-weighted inner product ( 6 b) follows when 
applying f_ 1 dxT m (x )... to both sides of ( 6 a). 


yields the original spectral function via A(tj) = 1 A(— + 
b ), where we omitted to specify the indices a, 6 , as in 
most of the rest of this paper. The Chebyshev moments 
for A(x) can be computed analogously to the Fourier 
coefficients (9) 

Pn = J dxA(x)T n (x) = {ifolifn), ( 12 a) 
IV’n) =T n (HMo). ( 12 b) 
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Inserting the recursive definition (4) of T n (fH) in the 
definition (12b) of \i/j n ), one obtains a practical calcula¬ 
tion scheme for the power series expansion of T n (fH), and 
by that for the Chebyshev states \ip n ) in (12): 

IV'n) = 2'H|Vv l —l) - \1pn-2), \ijJi} ='H\ijJo}- (13) 


C. Analytic continuation 

A comparison of (9b) and (12b) clarify why linear pre¬ 
diction is an equally justified approach for a Chebyshev 
as for a Fourier expansion. 

Rewriting the “evolution operators” that appear in 
(9b) and (12b) as 

exp>(—iriH a ,b=o ) and cos(n arccos("H a ,&)) (14) 

makes it clear that we deal with analytic functions of n if 
we consider n as a continuous complex variable. By (9a) 
and (12a), this makes the Fourier and the Chebyshev 
coefficients f(t n ) and p n analytic functions of n, too. 
In addition, Re exp (—iriH a ,b) = cos {nH a ,b) shows that 
the real part of the functional dependence of the time 
evolution operator on n is the same as for the “Cheby¬ 
shev evolution operator”, when neglecting a redefinition 
of oscillations frequencies. As this redefinition can be ac¬ 
counted for by the fitting procedure, the particular form 
of the surrogate function g{t) in (1) is equally suited to 
analytically continue both types of expansions. A fun¬ 
damental theorem from complex analysis then tells us 
that if linear prediction provides us with a function g(t) 
that locally agrees with f(t n ) or /i n , we know that this 
function globally agrees with f(t n ) or p n . Of course, in 
practice these arguments are to be taken with care, as we 
will never numerically find a function g(t n ) that agrees 
exactly with the local data {g n ,t n }. 


D. Outline of the paper 

We will first study the convergence properties of the 
Chebyshev expansion of discontinuous (spectral) func¬ 
tions in the thermodynamic limit. This allows to derive a 
new scheme for a Chebyshev series defintion that leads to 
exponential convergence and allows to reduce expansion 
orders in practical calculations by a factor ~ | (Sec. II). 
We then apply these results to the computation of spec¬ 
tral functions for finite systems (Sec. Ill), and discuss the 
relevance for matrix product state (MPS) based compu¬ 
tations (Sec. IV). After that, we describe the approximate 
equivalence of the Chebyshev recursion to time evolution 
and show how this leads to a novel time evolution algo¬ 
rithm (Sec. V). Finally, we conclude the paper (Sec. VI). 


II. SPECTRAL FUNCTIONS IN THE 
THERMODYNAMIC LIMIT 

A spectral function for a system of finite size L has 
a finite-size peak structure due to an agglomeration of 
eigenvalues that is not present in the thermodynamic 
limit. In a weakly interacting system, this agglomer¬ 
ation happens around the positions of the eigenvalues 
of the corresponding noninteracting (single-particle) sys¬ 
tem. This argument gives us the best, though still very 
rough, estimate W s i ng i e /L for the spacing of finite-size 
peaks, where Wsi ng i e is the single-particle bandwidth. At 
a much smaller spacing than that, spectral functions have 
an underlying delta peak structure, as is obvious from 
definition (8), which can be rewritten as 

A(u) = 6 ( w ~ ~ ^ef)), ( 15 ) 

i 

with weights Wi = \(ipo\Ei) | 2 . The delta peak structure 
merges to a (section-wise) smooth function only in the 
thermodynamic limit. 

Expanding the spectral function of a finite-size system 
in orthogonal polynomials is a very efficient way to not 
resolve either finite-size peaks or the delta peak structure, 
but to extract only the smooth function of the thermody¬ 
namic limit, as e.g. discussed in Ref. 1- . It is this function 
of the thermodynamic limit that we are interested in, and 
for which we start our discussion. 

A. Discontinuity of spectral functions 

The state |^o) and the energy E re f in (15) are gener¬ 
ally associated to the ground state of a certain symmetry 
sector N of H, which for fermions is typically a particle 
number. The reference energy for \if 0 ) = c*\Eq) then 
is the Fermi energy, which is the ground state energy 
E le { = Eq- 1 of the contiguous symmetry sector of \’ipo) 
(or E re f = Eq +1 for a hole excitation). The weights 
Wi = |(^o|^)| 2 in the spectral function (15) can be non¬ 
zero only for eigenstates | Ef) and eigenvalues Ei from the 
sector N. The particular meaning of E le f as a ground 
state energy then implies that even if the global spectral 
function 

AgiobaiM = ^2 S(uj - (Ei - E ref )) (16) 

i 

is smooth, the weights W t generally introduce a disconti¬ 
nuity at ui = 0 (we use the term global here, as j^o) usu¬ 
ally is a local excitation associated with a certain quan¬ 
tum number). 

B. Convergence of Chebyshev series expansions 

The convergence of the Chebyshev moments g, n —> 0 
of a function f{x) in the limit n —> 00 can be charac¬ 
terized by the degree of differentiability of f(x), similar 
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to a Fourier expansion. Let k denote the highest inte¬ 
ger for which the /cth derivative of f(x) is integrable: If 
/( x) is smooth (k = oo), the envelope of /z n converges 
exponentially to zero with respect to n; if f(x) is a step 
function (k = 1), the envelope converges algebraically 
with and if f[x) is a delta function, the envelope re¬ 
mains constant. In general, the order of convergence is 
at least Although in Ref. 15, this is stated for mo¬ 
ments computed with the weighted inner product (5b), 
it also holds for moments computed using (6b) (see Ap¬ 
pendix A). In practice, we are not interested in the limit 
n — > oo, but rather in intermediate values of n: but also 
here, the degree of differentiability of A(w) helps us to 
learn something about the convergence of /z„. 

Consider a typical discontinuous spectral function 
A > (uj) as shown in the top panel of Fig. 1. Its cor¬ 
responding Chebyshev moments /z> are computed by 
numerically integrating (6b) and shown in the bottom 
panel of Fig. 1 as blue circles. The blue line in the in¬ 
set shows the envelope of /z>, which evidently decreases 
algebraically to zero. 

Now note that continuity of A > (w) at w = 0 can easily 
be restored by defining 

A > (u) = A > (u)-A > ( 0). (17) 

The green crosses (lines) in the bottom panel of Fig. 1 
show that the Chebyshev moments /z> of A > (uj ) converge 
exponentially for the values of n considered in the plot, 
i.e., qualitatively differently than This is observed al¬ 
though A > (uj ) is not smooth, but only once differentiable 
(kink in first derivative at to = 0). 

While the construction of A > (u}) is completely general, 
for the particular case of a fermionic spectral function, 
another way of constructing a continuous function from 
A > (uj) has been favored: In the appendix of Ref. 17, it 
was mentioned that the Chebyshev expansion of the full 
spectral function 

A(uj) = A > (ijo) + A < {— w), (18) 

A*(u>) = {^\5(u - (H - E 0 )M), 

obtained by summing over particle (>) and hole (<) con¬ 
tributions, is much better suited for a Chebyshev expan¬ 
sion than A^(cc), as it lacks the discontinuity. In Ref. 13 
it was then pointed out that the full A(ui) is smooth and 
therefore, Chebyshev moments should decrease exponen¬ 
tially, which would allow to use linear prediction. In 
general, it is not true that A(oz) is smooth, due to the 
possibility of van Hove singularities, as appear e.g. for the 
U = 0 case of the spectral function of the single impu¬ 
rity Anderson model (SIAM) (see Appendix B, Fig. 11). 
Still, Afjj) is likely to be smooth, and for the present ex¬ 
ample, it is. The bottom panel of Fig. 1 therefore shows 
that the Chebyshev moments for A(w) decrease at the 
same exponential rate as the moments /z> of A > (w). 

The statements about the qualitatively different con¬ 
vergence behaviors of A > (uj), A > (uj) and A(az) are con¬ 
firmed for further typical examples in Appendix B. 



n/a (1/v) 

FIG. 1. (Color online) Top: Typical example of a discon¬ 
tinuous spectral function due to the restriction to a given 
symmetry sector. In this case, this is the particle contribu¬ 
tion of the spectral function of the single-impurity Anderson 
model (SIAM) with semi-elliptic bath density of states of half¬ 
bandwidth 2v and interaction U/v = 4, 1,1 taken from Ref. 

The spectral function is given by (8) using |^o) = c^Ao), 
where \Eo) denotes the half-filled ground state and E le f is the 
Fermi energy Eo. Here, only the shape of the scalar function 
is of importance, therefore we postpone the model definition 
to (30). The same spectral function is obtained for the local 
density of states of the first site for spinless fermions hop¬ 
ping on a semi-infinite chain with tunneling v and and an 
interaction of U/v = 4 that acts only on the first site. Bot¬ 
tom: Comparison of convergence of the Chebyshev moments 
of A y (uj) with its redefinitions A > ( w) and A(w), giving rise 
to moments gif, Jif and /z n , respectively. The full spectral 
function (18) for this example is A(ix) = A > ( ui) + A > (— to) as 
here, A < (uj) = A > (u>) due to particle-hole symmetry. All of 
this is for the setup b = 0 using a rescaling of a = 100u in 
(10). 


C. Comparison of setups b = 0 and b ~ —1 

In Ref. 14, we pointed out that the choice b = 0 in 
(10) is computationally much less efficient than the choice 
b ~ — 1 (called u b ~ —a” in Ref. 14). Whereas construct¬ 
ing a Chebyshev expansion of the full spectral function 
A(lo) requires choosing 6 = 0, this is not the case for 
A(co). For A{uj), we can therefore use the exponential 
rate of convergence to quantify the amount of spectral 
information that the Chebyshev recursion extracts from 
H in the setups 6 = 0 and 6 ~ —1, and by that under¬ 
stand the observations of Ref. 14 quantitatively. 
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FIG. 2. (Color online) Top: Integrand for computation of 
moments for the spectral function shown in Fig. 1. in the 
two setups 6 = 0 (left) and b = —0.995 (right). Bottom: 
Comparison of convergence of moments computed with the 
integrands shown in the top panels. In all of that, a = lOOu. 


The key observation to make is that the integral 

Vn = J i dx A> b (x)T n (x) (19) 

extracts a highly different amount of information about 
the structure of A > (co) depending on how a and b in (10) 
are chosen when generating A> b (x). 

Throughout the whole paper, we keep a = 100v fixed 
to guarantee the numerical stability of the Chebyshev re¬ 
cursion for the typical system sizes of around L > 80 that 
are large enough to display “thermodynamic limit behav¬ 
ior” . If we chose a smaller, we could only stably compute 
“small” systems or we would have to resort to the tech¬ 
nique of energy truncation , which is strongly prone to 
errors. Furthermore, in the MPS context, it is impor¬ 
tant to compare only computations in which a is kept 
constant: constant a means constant effective hopping 
energies - in and by that a constant amount of entan¬ 
glement production in a single iteration step of (13). The 
parameter 6, by contrast, can be chosen freely without af¬ 
fecting the numerical stability, and in principle, without 
affecting entanglement production in MPS computations. 

The top panels of Fig. 2 show the convolution of 
Aab ( x ) with Chebyshev polynomials T n {x) of different 
degree n for the two setups 6 = 0 and 6 = —0.995 ~ —1. 
The highly increased oscillation frequency that is evident 


in the setup 6 = —0.995 can be understood by looking at 
the natural stretching of the frequency scale of Cheby¬ 
shev polynomials close to the boundaries of [—1,1]. Ex¬ 
pressing the integral (19) by substituting x = cos 9 

r° 

//>=—/ d9A^ b (cos9) cos(n9) sin0, (20) 


one arrives at a convolution with the regularly oscillat¬ 
ing cos (nO). Consider now the interval of width 0.05 on 
[—1,1], which corresponds to the (single-particle) sup¬ 
port of A> b (x) in the example of Fig. 2. By comput¬ 
ing the integral widths under the map x = cos 0, one 
learns that placing the support in the “boundary region” 
[—0.995,-0.95], as results for 6 = —0.995, increases the 
resolution by a factor ~ 6.4 compared to placing it in 
the “center region” [0,0.05], as results for 6 = 0. These 
effects are well known boundary effects of the Chebyshev 
polynomials that are exploited also in the solution of dif¬ 
ferential equations. 

The bottom panel of Fig. 2 shows the Chebyshev mo¬ 
ments obtained in the 6 = —0.995 setup, for A > (cj) (blue 
circles) and for A > (ui) (green pluses). As mentioned be¬ 
fore, for this setup, no Chebyshev expansion of the full 
spectral function A(oS) is possible. Instead, we compare 
the 6 = —0.995 results to the 6 = 0 results, depicted as 
red crosses. It is evident that the Chebyshev expansion 
in the 6 = —0.995 setup converges much faster than the 
one in the 6 = 0 setup. After n = 1000 iterations, the 
magnitude differs by more than 100. 

This difference directly appears in the error of the 
Chebyshev series, as stated by the following general rule: 
the order of the error e of a Chebyshev (or Fourier) series 
representation of a function that is truncated at n = IV 
can be estimated by (see Ref. 15, Chap. 2.12) 

£ = O(iin) if Mn converges exponentially, (21a) 

£ = 0{N p, N ) if n n converges algebraically. (21b) 


D. Linear prediction for the Chebyshev expansion 


The main motivation for studying the convergence of 
different Chebyshev expansions in the previous sections 
lies in the possibility to extrapolate exponentially de¬ 
creasing sequences with linear prediction. As discussed in 
the introduction, the latter allows an extremely high gain 
in resolution, if its application is justified. For details on 
linear prediction, see Appendix C. 

In what follows, we compare the known approach of 
using linear prediction for the Chebyshev expansion of 
A(co), as suggested in Ref. 13, with the new approach of 
extrapolating the Chebyshev expansion of A> (w). 

We first compute the Chebyshev moments of the step 
function that has the discontinuity of A > (ui) at u = 0, 
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which transforms to x = — 6 for A > (x), as 


AC P = [ dxT n (x) 
J-b 


( 22 ) 


1 /cos[(n + 1) arccosx] 


;( 


n + 1 

cos[(n — 1) arccos x\' 1 


n — 1 


-b 


The Chebyshev moments of A > (x) are then given by 

Jl>=p>-A>( OK- (23) 


and are accessible by linear prediction, as they decrease 
exponentially. 

The core problem in this new approach is that the 
value A > (0) of the spectral function is, in general, un¬ 
known prior to linear prediction. But it fulfills the fol¬ 
lowing self-consistency problem, which can be iteratively 
solved: Choosing a start value Aq (0) for A> (0), we com¬ 
pute jH n , extrapolate the sequence up to convergence, and 
then use the extrapolated sequence to reconstruct A > (w), 
which provides us with a new value A> (0). We repeat 
the procedure until the new and the old version Af' (0) 
and Af, ^O) agree. This procedure is found to converge 
stably and quickly for all examples studied (see also Ap¬ 
pendix B). 

Fig. 3 compares the approach of reconstructing the full 
spectral funciton A(oj) from A > (u>), using linear predic¬ 
tion for the expansion of A(u) in the 6 = 0 setup, with 
the new approach of using linear prediction of A(uf) in the 
6 = —0.995 setup. We take the function of Fig. 1 as input 
function that shall be reconstructed. In the top panels 
of Fig. 3, we compare both setups for N = 200 com¬ 
puted moments that are then extrapolated to N 1000 
until they converge to a value of 10 -6 . We choose this 
comparatively small number of computed moments, as 
in MPS algorithms the number of moments that can be 
computed in a controlled way is strongly limited. 

The upper left panel of Fig. 3 shows that already for 
N = 200, our approach (dashed red line) allows a very 
good reconstruction of the input function. In the up¬ 
per right panel we show the error of this reconstruction, 
which becomes maximal at the second peak of the input 
function and is of order 10 —2 , i.e., a relative error of a 
few per cent. The situation is very different for the ex¬ 
trapolation scheme of the full A(ui) that uses the 6 = 0 
setup. For N = 200 computed moments, large errors are 
observed in both top panels of Fig. 3. 

In the bottom panel of Fig. 3, we plot the maximal 
error, defined as max^o |A reconst (u;) - A inpu t(w)|, ver¬ 
sus different values of the number of computed moments 
N. An orders of magnitude reduction of the error is seen 
upon using our over the previous approach. If one com¬ 
pares the expansion order N for which an error of 5- 10 -3 
is reached (N ~ 250 in the 6 ~ — 1 setup, and N = 1200 
in the 6 = 0 setup), one recovers the factor ~ 6 that has 
been derived in the previous section. 



N 

FIG. 3. (Color online) Top: Input spectral function and re¬ 
constructed spectral functions using linear prediction for N = 
200 computed Chebyshev moments. We compare our pro¬ 
posal with the original proposal, 1 where for the 6 = 0 setup 
the expansion of the full spectral function was extrapolated. 
Left: In the first case, we use the Chebyshev expansion of 
A > (u}) in the 6 = —0.995 setup (red dashed lines), and in the 
second, we use the Chebyshev expansion of the full spectral 
function A(to) in the 6 = 0 setup (solid green line). Right: 
Error of these functions |A reconst (u;) — Ai nput (cj) | for both se¬ 
tups. Bottom: The error max„> 0 | A reC onst(ta) — Ai nput (a;)| 
versus number N of computed Chebyshev moments. Here, we 
also show results for using the Chebyshev expansion of A > ( uj ) 
in the 6 = 0 setup (light green crossed line). This is differ¬ 
ent from using the full spectral function A(ui) in the 6 = 0 
setup. Lower error levels only occur for much higher higher 
expansion orders than shown in the panel. 


Only at very high expansion orders, that in practice 
can often not be reached by MPS computations, the orig¬ 
inal approach allows to reach smaller error levels for the 
presently studied generic example. More examples are 
studied in Appendix B. 


III. SPECTRAL FUNCTIONS FOR FINITE 
SYSTEMS 


Let us now study the case of finite systems, where a 
discretized representation of the spectral function is used 
for reconstruction. The general previous arguments are 
still valid, but several technical details have to be taken 
into account. In particular, we suggest a new discretiza¬ 
tion scheme suited for reconstruction with Chebyshev 
expansions. Such a discretization scheme can be used 
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for problems that allow to manipulate the discretization 
of the spectral function. This is e.g. the case for im¬ 
purity models, for which the discretization of the input 
bath spectral function determines the discretization of 
the spectral function. Still, the following discussion is 
also relevant to e.g. finite lattice models for which the 
discretization is physically constrained. 

To construct a discrete representation of a continu¬ 
ous function A(iS), we employ the scheme that is used 
to discretize the hybridization function of impurity mod¬ 
els in the numerical renormalization group. This pro¬ 
ceeds as follows. For L given discretization intervals 
[wz,wj+i], l = 1,...,L, we compute discrete weights V 2 
and eigenvalue positions e; by 

r^i+i 

V?= cLjA(uj), (24) 

Juji 

pU)l +1 

€l = 772 duiujA(u]). 

V l Jui 

The first line associates a weight and the second line 
a representative energy with an interval of energies 
[wzjWz+i]. For the energy, one could e.g. take the simple 
average |(wz + uii+i)- Eq. (24), by contrast, produces 
an average using the weighting function Aj-A(o;), which 

attributes more weight to peaks of A(w). 

We choose the left boundary of the first interval u>i and 
the right boundary of the last interval such that the 
distance ojl+i — wi is minimized but 

puiL+i r°° 

/ du A{uf) > 0.999 / duA(u), (25) 

J Uli J — oo 

where the integrand is non-negative, which guarantees 
that [wijWl+i] contains almost the complete support of 
A > (w), but minimizes finite-size effects. The intermedi¬ 
ate values of the discretization intervals { 0 J 2 , ...,wz,} can 
be chosen using a logarithmic discretization, as done in 
NRG. 8 But if an unbiased resolution is wanted, one usu¬ 
ally chooses a linear discretization 1 ’ 

coi = (l — l)Z\c 0 T cci, l = 2,..., L 

Auj = j(wl+i - wi). (26) 

As Chebyslrev polynomials do not show an unbiased 
energy resolution as they oscillate much quicker at the 
boundaries of [-1,1] than in the center, the linear dis¬ 
cretization will first resolve the finite size (discrete) struc¬ 
ture close to the boundaries of [—1,1]. We suggest to 
adapt the discretization to account for the cosine map¬ 
ping (20) of the energy scale that is responsible for this 
phenomenon. 

Let us study the case of even L (for odd L, see Ap¬ 
pendix D) and assume without loss of generality that we 
want as many intervals {wz,wz_|_i} on the positive half¬ 
axis as on the negative half-axis, which implies 

(27) 


As we already know uil+i from (25), we only have to 
fix the intermediate interval boundaries {u)l/ 2 + 1 , 

We define 

wl/ 2 +! = a(cos{6 L / 2 +lA0)-b), l = l,...,L/2 

A 9 =%(O L+ 1 - 0 L/2 ), 

Ql /2 = arccos b, 

9l +1 = arccos(6 + < ^ ± )- (28) 

Using these definitions, a discrete representation 
Adiscr(w) of A(oj) (in the sense that Adi scr (w) —> A(cj) 
for L —► 00 ) is given by 

AdiscrM = (lpo\S(u - H)\lpo), 

Hw = eiSw, l, l' = 1,..., L 
mi = Vi, (29) 

where H £ M. LxL and |^o) £ R L , and the parameters ei 
and Vi are given in (24). This is consistent with defini¬ 
tion (8) if we realize that this is a single-particle Hamil¬ 
tonian for a particle that is in either of e; energy states 
with probability V 2 . The reference energy would be the 
ground state energy of the vacuum E le f = 0. To obtain 
the step function behavior of A > (w), we project out the 
positive energy contributions from the initial state |^o)- 

In Fig. 4, we show the reconstruction of spectral func¬ 
tions based on the linear prediction of the moments com¬ 
puted for A^ iscl .(w) using the operator-valued Chebyshev 
expansion presented in Sec. IB for the “Hamiltonian” de¬ 
fined in (29). This is analogous to the top left panel of 
Fig. 3, which treated the thermodynamic limit. 

For the finite-size system, the specific choice of dis¬ 
cretization is important and we compare the linear and 
the cosine discretization in the top panels of Fig. 4 for 
the expansion order N = 200 and a system size L = 80. 
From the large error at w = 0 for the linear discretiza¬ 
tion (red dashed line) seen in the right top panel of Fig. 4, 
which was not present in the thermodynamic limit (red 
dashed line in right top panel of Fig. 3), we conclude 
that the linear discretization starts resolving finite-size 
features close to w = 0 already for N = 200. The lower 
panels then show how the error behaves as a function of 
the number of computed moments for different system 
sizes L. While the cosine discretization follows the error 
of the thermodynamic limit quite closely for low values 
of N and lattice sizes of L > 80, it almost saturates in 
a plateau for higher expansion orders, and only starts 
increasing slightly for very high expansion orders. For 
the linear discretization, neither the close correspondence 
with the thermodynamic limit is observed, nor does the 
error only moderately depend on the expansion order: 
Instead, the error increases exponentially for high values 
of N , as then, finite-size features are inhomogeneously 
resolved. Both features make it difficult to determine 
the value of N for which the computation of Chebyshev 
moments should be stopped in order to obtain a minimal 


w l /2 — 0 - 


error. 






N N 

FIG. 4. (Color online) Reconstruction of the spectral function 
of Fig. 1 represented by the discrete Hamiltonian (29). Top 
left: Input function and reconstructed functions using L = 
80, N = 200, a = lOOu, b = —0.995. We compare the linear 
discretization (26) with the cosine (28) discretization. Top 
right: Difference of input and reconstructed functions of the 
top left panel. Bottom left: The error max„>o |A reconst (u;) — 
j4in P ut(w)| versus number N of computed Chebyshev moments 
using a cosine discretization. Bottom right: Error for linear 
discretization. 


IV. IMPLICATIONS FOR MPS 
REPRESENTATIONS 

What is the relevance of the previous results for matrix 
product state (MPS) based computations of the spec¬ 
tral function for a given matrix product operator (MPO) 
HI Repeated MPO operations on MPS create entangle¬ 
ment, which eventually makes manipulating and storing 
MPS computationally very costly. Manipulations, such 
as applying H to states \t n ) in the recursion (13), or 
performing subsequent time evolution steps e~ lHAt , can 
therefore only be carried out up to a certain recursion 
order n or time t, before hitting an exponential wall in 
computation cost. For time evolution algorithms, this 
has long been known, - ’ but this also limits computa¬ 
tions using the Chebyshev recursion. 

In the following, we show that the method intro¬ 
duced in the previous sections outperforms the previous 
approach: it extracts more spectral information from 

H when creating the same amount of entanglement or, 
which is equivalent up to technical details of the algo¬ 
rithm, using the same computation time. 


As an example, we compute the spectral function of the 
single impurity Anderson model (SIAM), which serves 
as a common benchmark ■ ■ ’ ’ and is highly rele¬ 
vant as it is at the core of dynamical mean-field theory 
(DMFT). 24 - 27 

The Hamiltonian of the SIAM is given as, 

tfSiAM = H imp + H hath + H hyh , (30) 

-^imp = U (not 2) (^04 2) ’ 

L b 

-^bath = EE GC l c la, 

1=1 <y 
L b 

Hhyb = EE( Vicl a ci a + H.c.) . 

1=1 <j 

By a unitary transform effected by Lanczos tridiagonal- 
ization, this can be mapped on the so-called chain geom¬ 
etry. But as this leads to higher entanglement, we simply 
order bath states by their potential energy which directly 
gives a one-dimensional array that can be treated with 
MPS. ' We solve the model for the semi-elliptic bath 
density of states 

- ^ im/i( “ ) =’ <3i) 

which is discretized according to the procedure discussed 
in Sec. Ill, and then yields the parameters e; and V). It 
is important to realize that here, we discretize the bath 
hybridization function whereas in Sec. Ill, we discretized 
the spectral function. While Sec. Ill did this to illustrate 
the effect of discretization for a toy model for which the 
spectral function was known from the beginning, in the 
present case, a true many-body compuation is involved. 
In the present case, the relevant discretization parameter 

is the bath size Lb = L — 1, and no longer the system 

14 

size. 

We compute the spectral function (18) of the impu¬ 
rity Green’s function, where the initial states are single¬ 
particle excitations of the ground state: \ifo a ) = c o<rl-®o) 
and | ifoa) = co a \E 0 ). As we consider the particle-hole 
and spin-symmetric case of (30), we only need to com¬ 
pute one Chebyshev recursion; to be precise: |i/’o) = 
cj t | A 0 ). We compare our results with dynamic DMRG 
results from Ref. 16, which are believed to be highly re¬ 
liable. In particular, we compare computations in the 
formerly suggested setup 1 that uses the Chebyshev 
recursion for b = 0 in (10) and reconstructs the full spec¬ 
tral function A(uj ) using linear prediction, and the one 
suggested here, that uses b = —0.995 and reconstructs 
the shifted spectral function A(uj ) using linear prediction. 

In the top left panel of Fig. 5, we show computations of 
the spectral function of the SIAM for L = 80 for N = 260 
in the b = —0.995 setup, and N = 900 in the 6 = 0 setup 
and compare it with the result of Ref. 16. We choose 
these two expansion orders, as they lead to a comparable 
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FIG. 5. (Color online) Comparison of MPS computed spec¬ 
tral functions in the two setups studied in the previous sec¬ 
tions, with data by Raas et al. . Solid lines refer to the new 
method that uses A(w), dashed lines refer to the method that 
uses the full spectral function A(ui). Top left: For L = 80 
and two exemplary expansion orders. Top right: Errors of 
comparison in top left panel. Bottom left: Plot of the max¬ 
imum error versus expansion order N. Bottom right: Plot 
of the maximum error versus computer time. 


maximum error, as shown in the top right panel of Fig. 5. 
In the 5 = —0.995 setup, this maximum error is slightly 
smaller. Around to = 0, by contrast, the error in the 
b = —0.995 setup is much smaller. If we compare the 
computation time that is needed to reach this precision 
(max|A reconst (w) — Ai nput | ~ 0.015/v), we find that the 
b = —0.995 setup required ~ 145 min whereas the 5 = 0 
setup required ~434min. If one makes this comparison 
for a slightly larger error (max|A reconst (w) - A input | ~ 
0.025/u), realized for expansion order N = 120 for the 
b = —0.995 setup, and for expansion order N = 200 for 
the 5 = 0 setup, the comparison in computation times 
reads ~12min versus ~160min. 

When studying the convergence of the maximum er¬ 
ror with respect to expansion order N in the lower left 
panel of Fig. 5, we see that this is, after a sharp de¬ 
crease for low expansion orders, not monotonously de¬ 
creasing. The previously mentioned choices, N = 260 
in the 5 = —0.995 setup and N = 900 in the 5 = 0 
setup, both correspond to a minimum in the oscillations, 
as seen when inspecting the green solid (dashed) lines 
for the 5 = —0.995 (5 = 0) setup. The non-monotonicity 


makes general comparisons for the speedup difficult. But 
the lower right panel of Fig. 5 still shows that with only 
a few exceptions, the solid (5 = —0.995) lines are always 
clearly below the dashed (5 = 0) lines. The logarithmic 
abscissa therefore indicates a high speedup. 


V. COMPARISON TO TIME EVOLUTION 

It is interesting to compare the efficiency of the 
available MPS algorithms to extract spectral informa¬ 
tion from H. Candidates are, aside from the dy¬ 
namic DMRG, which is believed to be computationally 
highly costly, time evolution and recursive algorithms. 
The latter are, in particular, expansions in Chebyshev 
polynomials and the Lanczos algorithm. ’ Lanczos is 
numerical instable as the basis that it spans looses its or¬ 
thogonality for high numbers of iterations. This seems 
to to disqualify Lanczos as a high-performing candidate. 
Therefore the main question is whether the Chebyshev 
recursion can more efficiently extract spectral informa¬ 
tion from H than time evolution algorithms. 

To answer this question, in the following, we exploit 
the fact that for 5 = 0 in the limit a —> oo, the Cheby¬ 
shev expansion becomes a Fourier expansion, and the 
Chebyshev states directly describe the time evolved sys¬ 
tem. This is different from the procedure of computing 
the time dependence of a Green’s function via a Fourier 
transform of its spectral function. • “ For comparison, 
we summarize the latter technique in Appendix E. 

The approximate equivalence of the Chebyshev recur¬ 
sion to time evolution can be used as a novel time evo¬ 
lution algorithm. This is interesting as for long-range 
interacting Hamiltonians H, the MPO representation of 
e~ lHt is not available, or only approximately.' Although 
it is possible to use so-called Krylov algorithms for such 
problems, this requires some programming effort, and is 
in general believed to be numerically rather inefficient 
as compared ot other time evolution algorithms. Long- 
range interacting problems appear e.g. if mapping a two- 
dimensional system on a one-dimensional chain, or in the 
solution of a SIAM using a star geometry as in (30). 


A. Statement of approximate equivalence 

The time evolution of a state |/>o) 

I = exp(-iHt)\ip 0 ), (32) 

can be approximately linked to the sequence of Cheby¬ 
shev vectors generated by starting from j^o) as follows. 

Choose a reference energy E le f in (10) that is charac¬ 
teristic for the initial state of the time evolution and the 
Chebyshev recursion. When computing the time evo¬ 
lution of the Green’s function with \ipo) = c'\E 0 ), one 
chooses E le { = E 0 , if |/ 0 ) is not an eigenstate, we choose 

Ere f = {ipo\H[ijj 0 }. 
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Then define | <j>(t)) = exp(iE le {t)\ij)(t)) and H = (H — 
E re f)/a as in (10) in the 6 = 0 setup. Here a has the 
meaning of an inverse time step of unit energy. With 
these definitions, (32) reads as 

| </>(£)) = exp(—iaHt)\ipo) 

= (cos (aHt) — isin(a'Ht))\'ipo) 

= I0COS (t)) ~ #sin (t)). (33) 

Let us discretize time by defining f n = (C then 

l^cos(in)) = cos(nH)l^o), (34a) 

l^sin(in)) = sm(nH)\i>o). (34b) 

We now want to compute the action of cos(n'H)lV’o) on 
\ip 0 ) using a recursion that only involves the action of TL. 
This is not possible with the standard recursion for the 
cosine function, as shown in Appendix F 1. 

Let us instead consider the action of the Chebyshev 
polynomials 

T n (TL) = cos(narccos("H)) (35) 

on \ipo)- This action approximately reproduces the ac¬ 
tion of the plane cosine function, if we consider every 4th 
iteration, i.e. introduce the new index n' = 4n, n £ N: 

T n ' ('H)IV’o) = cos (n '(^ - 'H^\tpo) + e^OIV’o) (36) 

= cos (n'Hj |V’o) + e(n)\i/jo)- n = 0,4,8,... 

In the first line, we used the Taylor expansion 
arccos("H) = f — T~L + \l-L 3 + • • •, that leads to the error 
function e(n') (Appendix F 5), and in the second line, we 
used n' | = 2i m, n £ N, which obviously drops out of the 
argument of the cosine (also see Appendix F2). 

The error e(n') is bounded by (F16) 

\e(t n ')\ = - t — if t < ferr 
^err 

te rr = ^ (37) 

( 7 ° 

where a is the spectral width of the initial state |^q) = 
|^o) around E lef , 

a= max \E k - E lei \, (38) 

\E k )e\ip 0 ) 

where “| E k ) £ IV’o)” refers to the decomposition of the 
initial state in eigenstates \E k ) of Ti 

|V>o) = 5> fc | E k ). (39) 

k 

The “spectral width” a is usually small compared to rea¬ 
sonably high values of the inverse time step a. If one 
is unsure of whether a was chosen large enough, one re¬ 
runs a calculation with a higher value of a and checks 
convergence. 


We can now compute the time evolution 

l^cos (tn')) = l^cheb ifn' )) 4“ e{t n ' ), 
l^cheb(in')>=^'(^)l^o), n' = 0,4,8,... (40) 

via the recursion (13) 

l^cheb(^n)) — 27^|0 c heb(^n—l)) |^cheb(ln— 2 ))^ 

|0cheb(ti)> = n\ih), n = 0,1,2,.... (41) 



t=n/a (1/v) t=n/a (1/v) 


FIG. 6. (Color online) Left: Time evolution of a particle cre¬ 
ated on the first site of a chain of length L = 100 with hopping 
v = 1 that is l^o) = Cg|vac). We compare the time evolu¬ 
tion of the Green’s function iG > (t) = e lEot (shown 
as blue crosses) with the Chebyshev moments = (ipofyn) 
(shown as green dots) obtained when computing the recursion 
\ipn) = 2H / a\tpn-i) — | tp n — 2). Hopping amplitudes vi are 
obtained from the discretization of the spectral function of 
Fig. 1. A qualitatively equivalent behavior is obtained for a 
chain with homogeneous hopping vi = v. Right: Long-time 
evolution of the Green’s function (blue crosses), and differ¬ 
ence of the Green’s function and the Chebyshev moments (red 
dots). The horizontal dashed line marks the the prefactor of 
the error estimate (37), which is computed as a 3 /a 2 = 8T0 -4 . 


B. Numerical examples 

1. Single-particle computation for SIAM 

Fig. 6 shows the numerically exact time evolution of a 
single particle created on the first site |^o) = Cq\Eo) °f 
a chain of 100 lattice sites. The spectral width therefore 
is a = 2v. The left panel of Fig. 6 plots the Cheby¬ 
shev moments /z n = (ifo\tn) obtained with a = 100c and 
the time evolution of the corresponding Green’s function 
iG > (t ) = e lEot {ifo\ V’(i)) for the time step 4. Every forth 
Chebyshev moment agrees with a value of the Green’s 
function. In the right panel, we show the long time be¬ 
havior of G > (t) and the difference G > (4n/a ) — /i± n . The 
difference is seen to be clearly below the conservative up¬ 
per bound (37), it remains of the order of a 3 /a 2 = 8-10~ 4 
up to very high times that correspond to 100 hopping 
processes (t = 100/c). 
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2. MPS computation for SIAM 

We now study the time evolution of the SIAM (7) 
in the star geometry for the single-particle excitation 
IV’o) = cl t \E 0 ) of the half-filled ground state \E 0 ). The 
left panel shows the time evolution of the corresponding 
Green’s function, computed with an MPS Krylov algo¬ 
rithm that imposes the error bound 

11 | i/j{t + At)) - exp(-iHAt)\ip(t)) \ \ < e kry , 

for a time step of At = A An error bound of £ kry = 
5 • 10 -4 suffices to reliably compute times up to 15/i>. 

In the MPS implementation of the Chebyshev recur¬ 
sion, we fix the global truncation error per iteration step, 
as discussed in Ref. 14 

| I | tn) — ( 2'H\t n -i ) — |tn— 2 )) I I < £ c he- (42) 

To achieve this, two options are available. If during the 
variational compression of (2"H|t n _i) — |f n _ 2 )), the trun¬ 
cation error exceeds £ c he, even when choosing a better 
and better guess state, one can either directly increase 
the bond dimension, or reduce the truncated weight per 
bond, which indirectly increases the bond dimension. 
While for the setup in Ref. 14, there were reasons to 
choose the former option, here we choose the latter as 
our Krylov algorithm uses a similar adaption. 

We compare the results of the Krylov algorithm with 
the Chebyshev algorithm (40). In the top left panel 
of Fig. 7, we plot the Green’s function iG > (t ) = 
e lE ° t (ipo\'il>(t)}. If imposing the same error tolerance 
£ c he = £kryi we obtain agreement of both algorithms only 
for short times. Only a much smaller tolerance for the 
Chebyshev algorithm £ c h e = j^£kry leads to agreement 
also for long times. We conclude that error accumula¬ 
tion in the Chebyshev recursion is much worse condi¬ 
tioned than in the time evolution algorithm, and even 
worse than what could be expected from the four “aux¬ 
iliary steps” made in (41) between each “physical time 
step”: imposing a tolerance £ c h e = j£kry for the Cheby¬ 
shev recursion is not sufficient to produce comparable 
results. 

The reduced error tolerance £ c he = pj^kry for the 
Chebyshev recursion comes at the price of an order of 
magnitude increase in the bond dimension compared to 
the Krylov algorithm, as shown in the top right panel 
of Fig. 7. But also for £ c h e = £kry, the Chebyshev re¬ 
cursion needs higher bond dimensions than the Krylov 
algorithm. The lower left panel compares the overlap of 
the Chebyshev-evolved and the Krylov-evolved states by 
plotting |1 - (^cheb|^kry)/(^cheb|V’cheb)|- With only few 
exceptions, this quantity is bounded by the theoretical 
prediction of (37), when setting a = 4v = U. The ex¬ 
ceptions are artifacts of the detailed implementation of 
the algorithms as the key observable G > (t ) is correctly 
computed, but still their existence suggest that the imple¬ 
mentation can be improved. Ignoring these exceptions, 
we see that the normalized overlap (V’chebl'^kry) deviates 



FIG. 7. (Color online) Time evolution of the single parti¬ 
cle excitation c^l-Eo) in the half-filled single impurity Ander¬ 
son Model (30) with semi-elliptic density of states of half¬ 
bandwidth 2v and interaction U/v = 4 for L = 40. Compu¬ 
tations using an MPS Krylov algorithm with error tolerance 
e kry = 5 ■ 10 -4 and the Chebyshev recursion (40) for different 
error tolerances £ c he = 5 • 10~ 4 and £ c h e = 5 • 10 -5 . Top 
left: Time evolution of Greens’s function. Both algorithms 
produce the same result upon using the smaller error toler¬ 
ance for the Chebyshev algorithm. The legend is found in 
the top right panel. Top right: Maximal bond dimension, 
located at the central bond. The Krylov time evolution leads 
to a much smaller maximal bond dimension, as its computa¬ 
tion produces a faithful result already with the relatively high 
error tolerance of £ kl . y = 5 ■ 10 —4 , for which the Chebyshev al¬ 
gorithm shows strong errors in the Green’s function. Bottom 
left: Difference of overlap of Chebyshev and Krylov evolved 
states, comparing the £ c h e = 5 • 10~ 5 with the £ kry = 5 ■ 10 -4 
computation. The difference of overlap is bounded by the 
analytical prediction of (37), except for few exceptions that 
lie above it. These exceptions are of purely numerical origin 
as they are not visible in any other quantity. For the high¬ 
est times shown, truncation errors have accumulated so much 
that the analytical prediction starts to fail. Bottom right: 
Bond entanglement entropy at center bond. 


from one only by a few percent even for long times. But 
these few percent come with a considerable growth of the 
entanglement entropy, as can be concluded by inspect¬ 
ing the lower right panel of Fig. 7. There, already the 
£ c h e = £kry case shows a considerably increased entropy. 

Aside from the two preceding fundamental reasons 
(different error accumulation, small difference of states), 
the increased bond dimensions in the Chebyshev algo- 
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rithm can also be related to a purely technical question: 
the variational compression in each Clrebyshev itera¬ 
tion produces a state that fulfills (42), but might be a 
state with unnecessarily high bond dimension m. Simi¬ 
larly to the DMRG ground state optimization algorithm, 
also variational compression can get stuck in local min¬ 
ima. Currently, we use White’s mixing factor to avoid 
this. A recent publication suggests an even better strat¬ 
egy and explains these problems concisely. 



gled basis of this subspace than a time evolution algo¬ 
rithm: it extracts less spectral information when fixing 
a maximal entanglement entropy. But these arguments 
directly hold only for the u b = 0 setup” of the Cheby- 
shev recursion, in which it is transparently comparable 
with a time evolution algorithm as there is a one-to-one 
correspondence of time evolution steps and iterations of 
the recursion. 

Sections II-IV of this paper showed that the b = 0 
setup is the computationally least favorable setup of the 
Chebyshev recursion, and a b ~ — 1 setup much better. 
Still the gains in computation time of the b ~ — 1 setup 
over the b = 0 setup shown in the bottom right panel 
of Fig. 5 seem not to be sufficient to compensate the 
clear inferiority of the Chebyshev method shown in the 
upper right panel of Fig. 7. A definitive statement is 
difficult due to the non-monotonic behavior of the error 
in Fig. 5 and due to the fact that such a comparison is 
strongly affected by the details of the implementation of 
the algorithms, and not only by the principle nature of 
how strongly entangled its resulting basis states are. For 
this reason, the discussion on the most efficient method 
for computing spectral functions using MPS cannot be 
generally considered settled. 



FIG. 8. (Color online) Time evolution of a one dimensional 
fermionic Hubbard model on L = 90 sites, with an interaction 
of U/v = 4 and nearest neighbor hopping v starting from a 
product state (double occupation in the center of the system). 
Left: Chebyshev computation using £ c h e = 0.0001. Right: 
Krylov computation using £k ry = £ c he- Bottom: Detailed 
comparison for the occupation of specific sites. Chebyshev 
results are shown as dashed lines, Krylov results are shown as 
solid lines. Deviations are smaller than 1%. 

In general the subspace of the Hilbert space that is 
needed to be faithfully described in order to measure the 
spectral function can be spanned using different basis 
states. In principle, the most efficient spanning would be 
provided by the Lanczos algorithm, as the latter provides 
orthogonal states. But it is impractical due to numerical 
instability. The basis states provided by time evolution 
or the Chebyshev recursion are not orthogonal to each 
other, but can be stably generated. The numerical evi¬ 
dence discussed in the previous paragraphs indicates that 
the Chebyshev recursion generates a much higher entan- 


3. Expansion in Hubbard model 

Finally, we study the time evolution of the one- 
dimensional Hubbard model 

F Hubbard = [/ ^( n;t _l ) ( „ 4 _1) 
l 

-«E( c l c i+ ia+h-c.), (43) 

la 

starting from a product state with doubly occupied sites 
in the center of the system, and evolving this state at 
interaction U/v = 4, as shown in Fig. 8. We obtain 
very good agreement of the Krylov and the Chebyshev 
algorithm, although there is no rigorous a priori reason, 
for which the initial product state should have a narrow 
spectral width, i.e. small a in the sense of (38), as was 
the case for the single particle excited initial state. On 
the other hand, for example in the many studies on the 
eigenstate thermalization hypothesis >s ’ it is a frequently 
met assumption that for typical initial states the energy 
distribution around its mean value is extremely narrow, 
with a width of the order of the single-particle energy 
scale (see e.g. Ref. 40 Fig. 3b). 

VI. CONCLUSION 

We started by linking linear prediction to analytic con¬ 
tinuation, which explains why linear prediction is a rea¬ 
sonable method to extrapolate both Fourier and Cheby¬ 
shev expansions of spectral functions. In order to apply 
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linear prediction, we introduced a new method to avoid 
the algebraic convergence of Chebyshev moments (expan¬ 
sion coefficients) of generic step-like spectral functions. 
This amounts to a particular redefinition of the series ex¬ 
pansion that is based on a subtraction of the Chebyshev 
moments of a self-consistently determined step function. 

We then showed that this allows to reduce the expan¬ 
sion order by a factor ^ as compared to the existing 
method. For linearly scaling algorithms, as in exact 
diagonalization, 1,11 this means a reduction of computa¬ 
tion time of the same factor. But also for matrix prod¬ 
uct state computations high speedups are obtained. Fur¬ 
thermore, we showed how to adapt the discretization of 
hybridization functions of impurity models to the Cheby¬ 
shev method. 

Finally, we showed the approximate equivalence of the 
Chebyshev recursion to time evolution in a certain limit. 
This lead to a novel time evolution algorithm and al¬ 
lowed to transparently compare standard time evolution 
and the Chebyshev recursion in how efficient they ex¬ 
tract spectral information from an operator H. For ex¬ 
act representations, the Chebyshev recursion is superior 
to time evolution as the latter is equivalent to the least 
favorable setup of the Chebyshev expansion, which can 
be improved by the previously mentioned factor |. For 
matrix product state representations, our results indicate 
that the Chebyshev expansion is inferior: we observe a 
much higher entanglement production in the Chebyshev 
recursion than in standard time evolution. We identify 
as main reason for this an unfavorable error accumula¬ 
tion in the Chebyshev recursion that requires computa¬ 
tions at higher accuracy. So while in the history of the 
solution of differential equations for non-time periodic 
problems, Chebyshev expansions replaced Fourier expan¬ 
sions in the course of time, 1 ' in the matrix product state 
context such a transition now seems unlikely. Still the 
Chebyshev recursion provides an easy-to-implement and 
straight-forward way to compute spectral functions. 


APPENDIX 


Appendix A: Convergence speed 


Analogously to Ref. 15, Chap. 2.9, we give the argu¬ 
ment for the speed of convergence of the Chebyshev se¬ 
quence, computed with the non-weighted inner product 
of (6b) 


Hn = J dxf(x)T n (x) 

cLOf (cos 9) cos (nO) sin 9 

= Re [ d9 f (cos 9) e in9 , (Al) 

Jo 

where f(9 ) = /(cos 9) sin 9. We can then do k partial 
integrations, if f(9) is k times differentiable 



9"n — Re 


t=i 


G 


d9f {k \9)e ine 


(A2) 


where f^(9) denotes the jth derivative of f(9). If 
0) = 7r) = 0 for j = 0,..., k — 1, which is ful¬ 

filled for typical single-particle spectral functions as in 
Fig. 2, and if f ( - k \9) is integrable, (A2) constitutes an 
upper bound 0 (for the sequence /*„• 


Appendix B: Examples for linear prediction of 
Chebyshev expansions 


Relevant applications of the results of this paper are 
the computation of conductivities, 4 the computation of 
time evolution of long-range interacting systems,and in 
particular, the challenging solution of dynamical mean- 
field theory. 1 1 1 For example, the latter can usually not 
be accessed by combining analytical and numeric tech¬ 
niques as recently done for the Hubbard model in Ref. 42. 
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In Sec. II, we compared the reconstruction of a spec¬ 
tral functions using its extrapolated (linearly predicted) 
Chebyshev expansion. We focussed on a typical example 
for this discussion, given by the U/v = 4 spectral function 
of the half-filled SIAM with semi-elliptic bath density of 
states, which is shown in the top panel of Fig. 1. 

In this appendix, we support the arguments of Sec. II 
by showing further generic examples. Starting from a 
step-like input function A > (cu), we again compare the two 
reconstructions based on (i) linearly predicting the “sub¬ 
tracted” spectral function A > (u >) of (17) and (ii) linearly 
predicting the “full” (summed particle and hole contri¬ 
butions) spectral function A(oj) of (18). 

To consider generic cases, we study functions that show 
“features” at oj = 0 and at some distance, of order of 
the single-particle bandwidth, away from it. The most 
natural choice for constructing such functions are super- 
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positions of (non-normalized) Lorentzians and Gaussians 


f> , , __ / 0 for w < 0 

f { £ Wo e{ 0 , 4 } M^ 0 ) else, 

= (w - J 0 ) 2 + ’ 

(-J —OJn ) 2 

h, g (uj,uj 0 )=e . 


The function / > (w) is plotted for both choices in the top 
panels of Figs. 9 and 10 for r] = 0.2. 

Based on the same argument as is the basis for Sec. V 
in this paper (approximate equivalence of Fourier and 
Chebyshev expansion), Ref. 14 showed the decrease of 
Chebyshev moments for superpositions of Lorentzians 
and Gaussians, to be approximately exponential and 
oc e~ an , respectively. This behavior is observed in both 
center panels of Figs. 9 and 10. For high values of n, in 
Fig. 10, the decrease oc e~ an transitions into an expo¬ 
nential decrease, which is not in contradiction with the 
result of Ref. 14. This, as the discussion of Sec. V, make 
statements for intermediate values of n. 

In the bottom panels of Figs. 9 and 10 we then show 
the error obtained for the different methods of recon¬ 
struction. Linearly predicting A > (oj) yields considerably 
lower errors than linearly predicting the “full” spectral 
function A(u>). Only in the case of Lorentzians (lower 
panel of Fig. 9), using A(w) leads to lower errors for val¬ 
ues of n. 

Finally, in Fig. 11, we show results for the spectral 
function of the half-filled non-interacting SIAM with 
semi-elliptic bath density of states, which itself is semi- 
elliptic, 




0 for io < 0, 

0 for w > 2v, 



(B2) 


A > (uj) has a kink at ui = 2v, as can be seen in the 
top panel of Fig. 11. Therefore, Chebyshev moments 
decrease only algebraically, as seen in the center panel of 
Fig. 11. 

If we consider the error of the linear-prediction-based 
reconstructed A > (w), shown in the bottom panel of 
Fig. 11, we see that this yields much better results than 
the estimate (21) gives for a plain truncation of an alge¬ 
braically decreasing series. Concerning the comparison 
between the two methods of reconstructing A > (oj) and 
the full A(u), we observe that A > { w) yields to smaller 
errors throughout. 

We finally note that studying the spectral function of 
the non-interacting SIAM with a constant bath density 
of states as in Ref. 17 and Ref. 13, does not constitute 
a more general case. The analytic expression for this is 
very close to a single Lorentzian. 




N 

FIG. 9. (Color online) Top: Test function consisting of two 
Lorentzians (Bl). Center: Corresponding Chebyshev mo¬ 
ments in the three different setups p > , ju > and p,, analogously 
to the bottom panel of Fig. 2. Bottom: Error of recon¬ 
structed spectral function, analogously to the bottom panel 
of Fig. 3. All of this is for a = 100. 


Appendix C: Linear prediction 

In the context of time evolution linear prediction has 
been long established in the DMRG community, 11,1 but 
it has only recently been applied to the computation of 
Chebyshev moments. 1,1 The optimization problem for 
the sequence p n becomes linear, if the sequence can be 
defined recursively 


p 

Pn ~ ^ ^ QiPn—ii (Cl) 

»=1 

which is easily found to be equivalent to (1). 1J The strat¬ 
egy is then as follows. Compute N c Chebyshev moments, 
and predict moments for higher values of n using (Cl). 

The coefficients cq are optimized by minimizing the 
least-square error X^eA/t-t n ~ Mnl 2 f° r a su bset A/g t = 
{N c — ns t ,..., iV c — 1, N c } of the computed data. We 
confirmed ng t = iV c /4 to be a robust choice, (i) small 
enough to go beyond complicated low-order (short-time) 
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FIG. 10. (Color online) Top: Test function consisting of 
two Gaussians (Bl). Center: Corresponding Chebyshev mo¬ 
ments in the three different setups gP, ju > and g, analogously 
to the bottom panel of Fig. 2. Bottom: Error of recon¬ 
structed spectral function, analogously to the bottom panel 
of Fig. 3. All of this is for a = 100. 


FIG. 11. (Color online) Top: Particle spectral function 
of half-filled non-interacting SIAM. Center: Corresponding 
Chebyshev moments in the three different setups g > , g > and 
g, analogously to the bottom panel of Fig. 2. Bottom: Error 
of reconstructed spectral function, analogously to the bottom 
panel of Fig. 3. 


behavior and (ii) large enough to have a good statistics 
for the fit. Earlier, we chose ng t = IV c /2, which leads 
to a better statistics for the fit. But this improvement is 
not important, as we do not deal with stochastic data. 
Minimization yields 


Ra = —r, a = —R 1 r, (C2) 

Rji = ^ ' [in—jUn—ii Tj = ^ ' Mri—j/Ai- 

neA/tit nEA/fit 


Linear prediction is more prone to overfitting if choosing 
p to be very high. Therefore one should restrict the num¬ 
ber of coefficients to p = min(nfit/2,100). Furthermore, 
one adds a small constant <5 = 10 -6 to the diagonal of R 
in order to enable the inversion of the singular matrix R. 


Defining 

/ — ai — a2 —a 3 ... — a p \ 

1 0 0 ... 0 


\0 0 ... 1 0 / 

one obtains the predicted moments fiN 0 +n = (M n HN a ), 
where pn 0 = (mat c -i Mat c -2 • ■ • t l N c - P ) T ■ The matrix 
M might have eigenvalues with absolute value larger than 
1, either due to numerical inaccuracies or due to the fact 
that linear prediction cannot be applied as p n rather in¬ 
creases than decreases on the training subset A/fit- In or¬ 
der to obtain a convergent prediction, we set the weights 
that correspond to these eigenvalues to zero measuring 
the ratio of the associated discarded weight compared to 
the total weight. If this ratio is higher than a few per¬ 
cent, we conclude that linear prediction cannot be ap¬ 
plied. One can then restart the computation to increase 
the number of computed moments N c , and try applying 
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linear prediction for a higher number of moments. 


on integrals 


Appendix D: Discretization for odd L 

In the case of odd L we cannot equate the central inter¬ 
val boundary with 0 as in (27). Instead we have to choose 
the correct width for the “central interval” by choosing 
the neighboring boundaries W(£ +1 )/ 2 and ^(l+\)/ 2 +i cor¬ 
rectly. This is achieved by subtracting an offset Aw from 
each of the positive boundaries defined in (28), such that 
for l = 1,2,... 


dxw n (x)T n (x)e la ^ x b 


(2 ffno)c 


— iabt 


dx 


T n { x)t 


—oo 7T 

— iabt 


Vi ~ ■ 


= (-i) n (2-5 n0 )e- tabt J n (at). 


(E3) 


Appendix F: Comparison to time evolution 
1. Standard recursion for cosine 


u} (L+i)/ 2 +i — a(cos(0( i+ i)/ 2 +i + lAO) — b — AlS) , 

AO = f(0L+l - 0(L+1)/ 2+l), 

0 (l+i)/ 2 + i = arccos b, 

0l+ i = arccos(6 + ^y~), (Dl) 

Aw = ^ (cos(0 (L+1)/2+ i + AO) - b). 

For negative boundaries Aw has to be added instead of 
subtracted. 


As usual for a vector space of orthogonal polynomials, 
the space of cosine functions {cos(nx)}, where n £ N, x £ 
R, can be generated using a three-term recursion formula 

cos (nx) = 2cos(x) cos((n — l)x) — cos((n — 2)x), (FI) 

which can be proven using addition theorems. 

Rewriting this in the operator-valued form for the ar¬ 
gument "H, and acting on l^o) yields for the definition of 
\4>co S (t n )} in Eq. (34) 


|0cos(£n)) = 2 COS T~L 1 0 C os (tra— 1) ) - |0cos(^n-2))- (F2) 

Appendix E: Time evolution by Fourier transform 

But this provides no solution for our problem as the 
action of cos Oi on |</> CO s(Ui-i)) is not known. 

Given the Chebyshev expansion in frequency space 


A > (w) = -A > (- + b), 
a v a 7 

A>(x) = ^2w n (x)n n T n (x), (El) 


we can obtain the time evolution of a single Green func¬ 
tion, but not for the whole system state, by Fourier trans¬ 
forming 


G>(t)= dwA>(w)e iut 
J — OO 
i r°° 

= -/ dwA > (- + b)e iut 
a J -oo a 

/ OO 

dx A> {x)e ia ^~ b)t 

-OO 

/ OO 

dx y w n (x)n n T n {x)e la ( x ~ b)t 

n 

/ OO 

dx w n {x)T n {x)e ia(x ~ b)t 

-oo 

= e~ iabt y(2 - S n0 )(-i) n ii n J n (at ), 


(E2) 


where the last step (interchanging sum and integral) is 
only possible if the sum is absolutely convergent, or finite. 
The Fourier transform can be looked up in a handbook 


2. Shifted cosine 


Using an addition theorem 
cos (n(f — a)) = 

= cos {n |) cos(na) + sin (re|) sin(na) 

cos (na) 
sin (na) 

(cos(na) + sin(na)) else. 


(F3) 


if f e N 
if £ N 


Eq. (36) follows if setting a = H. 


3. Sine term 

Changing the initial conditions of (F4a) generates the 
polynomials T' n = sin(narccos(x)) 

T' n {x) = 2xT' n _ 1 (x) - r n _ 2 (x) t (F4a) 

T[(x) = V 1-x 2 , T'(x) = 0, (F4b) 

which approximates the sine function, in the same way 
as (36) approximates the cosine 

T' n (H) ~ sin (W) if n/4 £ N. (F5) 
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4. Bound for Arccos 


time-evolved state. Therefore, 


Bounding the Arccos works as follows 

„2n+l 


arccos 




(tv 


n—0 


4 n (2n + 1) 


= - - x - r(x) 


'•(x)=x 3 (- + |^ 


oo ( 2 n )x2n - 2 


4 n (2n + 1) 


(F6) 


OO P™) 

Using arcsin(l) = 4n ( 2 Z+i) = f we can bound 


r x = x“ - 


G 


°o /2n\ 2n—2 


<|IH ,6 


l(s 


l(s 


£ 

n—2 

oo 

E 

n—2 


o 


4 n (2n + 1) 

( 2 ;) 


4 n (2n+ 1) 


"'■‘“dr 1 -., 


1 ( 5-0 


= |£ II - - 
2. 3 

< 3 |x 


(F7) 


5. Error computation 


^n'(wfc) = cos(n'arccos(wfe)) 

= cos(n'wjt — nVfc) as n' = 4n, see Eq. (F3) 
= cos (n'ojk) cos(nVfc) + sin(n , Wfc) sin(mrfc) 

= cos (n'wfc) + efc(n / ), 
efe(n / ) = cos(n / Wfe)(cos(wrfe) — 1) 

+ sin(?r / o;fc) sin(mrfe) 

Up to here everything was exact. 

Now we can strictly bound the absolute value of the 

2 

error term using |sin(rfc)| < |ffc| and | cos(rfc) — 1| < ^ 
and trivially bounding cos (n'u>k) and sin(n'wfc) by one, 

|e fe (n')l < \n' 2 rl +n’\r k \ (Fll) 

Let us now define the energy eigenvector |U max ) for which 
the error becomes maximal, which is the one for which 
Wmax is maximal, i.e. 

Wmax = max Wfc (F12) 

|£Ve|bo> 

with which we compute r max and e max . We can then 
simplify further 


The approximation in (36) is based on the Taylor ex¬ 
pansion 

a,Tccos{n) = ^-n + ln 3 + o(n 5 ), (fs) 
2 6 

. H - U ref 

ri , 

a 

which reflects the fact that the arcus cosine is well ap¬ 
proximated already by the leading linear term around 
x = 0. 

The approximation of \<j>(t n )) that has been generated 
in this way, is good if a is large enough and becomes exact 
for a V oo. But how large does one have to choose a in 
practice in order for e(t n i) to be bounded by the wished 
accuracy? 

Consider the decomposition of the initial state in eigen 
states | E n ) of T~L 

|V>o) = 5> fc |i? fe ), (F9) 

k 

and, defining W& = Ek ~ gref , therefore 

T n (HMo) = ^c fc T n ( Wfe ) \E k ). (F10) 

k 

We are now only interested in indices n' that are mul¬ 
tiples of 4 n, as only those have the interpretation of a 


^ ^ Cfc V(^ ) | Ek} £max(?7 ) ^ ^ C/c \E k } — £max(?7 ) l^o) ■ 
k k 

We therefore arrive at 

Tn'i'H)\ipo) = cos(n , 'H)|V'o) + e(xi') l^o) (F13) 

The value of w max is determined by the cutoff of the 
distribution of eigenvectors | E k ) in |t/> 0 )- This can be 
a strict cutoff or a few standard deviations of Gaussian 
distribution, beyond which no contributions with numer¬ 
ically measurable weight occur. Let denote this cutoff or 
width a and define it analogously to w max , i.e. 

cr = max | E k - E le{ \ V 

^max (F14) 

\E k )e\ip 0 ) a 

If, e.g., I^o) is constructed by applying a single-particle 
operator to an eigenstate (e.g. the ground state) of H, 
a is the single particle bandwidth Wsi ng ie times a small 
factor of order 1. 

Finally we need to bound the error term r k (Appendix 
F 4) 

fk = + 0(u>l), \r k \ < gl w fel- (F15) 
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Using the definition of a, let us now bound e(n') 

Kn')| < 

< 

< 


< 


Or expressing this in units of time 

\e{t n ')\ = 7 T t <i t err 

^err 

a 2 

terr = ~~ (F16) 

(J° 

Inserting typical values, where v is a hopping energy 
for a single particle process: a = 2v, a = lOOv one finds 
terr = The accumulated error e{t n ') therefore re¬ 

mains smaller than 10~ 2 if t < 12.5-. 

V 


c(nOI 


P' 2f, L + »'l r n 


^' 2 (|) 2 < ax 




3 ,2/cr 
—n - 
2 3 


© 2 (D 

© 3 


+ ©g|wLxl 

1 ,2/<t 

+ « »(- 
3 \ a 

if n f < n' 
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